Manuscript_CAPTAIN_analyses

Author

Théophile L. Mouton

Published

November 5, 2025

Load libraries

Code
library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)
library(readr)
library(tidyr)
library(ggpubr)
library(betareg)
library(broom)
library(margins)
library(patchwork)
library(dplyr)
library(sf)
library(egg)
library(spatstat)
library(sf)
library(tidyverse)      # For data manipulation and visualization
library(rnaturalearthdata) # Additional natural earth data
library(biscale)        # For bivariate mapping
library(gridExtra)      # For arranging multiple plots
library(grid)           # For text elements in plots
library(colorspace)     # For color manipulation

Load the data

Code
# Read the CAPTAIN2 EDGE2 RDS file
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 FUSE RDS file
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 IUCN RDS file
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
# Use the same raster file for consistency
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Read the protected range fractions RDS file
protected_fractions <- readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))

# Read the continental shark conservation metrics CSV file
shark_metrics <- read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))

# Load fishing data
load(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))

EDGE2 continental 0.1 budget

Code
# Analyze non-zero cells in the prioritization output
cat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")
Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:
Code
total_cells <- nrow(CAPTAIN2_EDGE2_data)
nonzero_cells <- sum(CAPTAIN2_EDGE2_data$Priority > 0, na.rm = TRUE)
zero_cells <- sum(CAPTAIN2_EDGE2_data$Priority == 0, na.rm = TRUE)
na_cells <- sum(is.na(CAPTAIN2_EDGE2_data$Priority))

cat("Total cells in grid:", total_cells, "\n")
Total cells in grid: 232560 
Code
cat("Cells with non-zero priority:", nonzero_cells, " (", round(nonzero_cells/total_cells*100, 2), "%)\n", sep="")
Cells with non-zero priority:5618 (2.42%)
Code
cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")
Cells with zero priority:226942 (97.58%)
Code
cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")
Cells with NA priority:0 (0%)
Code
# Summary statistics of priority values
priority_summary <- summary(CAPTAIN2_EDGE2_data$Priority)
cat("\nSummary statistics of priority values:\n")

Summary statistics of priority values:
Code
print(priority_summary)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00000 0.00000 0.02331 0.00000 1.00000 
Code
# Distribution of non-zero priority values
nonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority > 0]
cat("\nDistribution of non-zero priority values:\n")

Distribution of non-zero priority values:
Code
priority_quantiles <- quantile(nonzero_priority, probs = seq(0, 1, 0.1), na.rm = TRUE)
print(priority_quantiles)
  0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
0.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 
Code
# Load one of your input raster files to extract the correct grid structure
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 data based on PUID
CAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_EDGE2_sf <- st_as_sf(
  CAPTAIN2_EDGE2_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_EDGE2 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_EDGE2 <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_EDGE2 <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_EDGE2_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the plot
print(CAPTAIN2_EDGE2_plot)

Code
# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),
  plot = CAPTAIN2_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

FUSE continental 0.1 budget

Code
# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 FUSE data based on PUID
CAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_FUSE_sf <- st_as_sf(
  CAPTAIN2_FUSE_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_FUSE <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_FUSE <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_FUSE <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank())      # Add this line

# Create the plot
CAPTAIN2_FUSE_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE

# Display the plot
print(CAPTAIN2_FUSE_plot)

Code
# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),
  plot = CAPTAIN2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

IUCN continental 0.1 budget

Code
# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 IUCN data based on PUID
CAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_IUCN_sf <- st_as_sf(
  CAPTAIN2_IUCN_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_IUCN <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_IUCN <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_IUCN <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_IUCN_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN

# Display the plot
print(CAPTAIN2_IUCN_plot)

Code
# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),
  plot = CAPTAIN2_IUCN_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Difference maps

Code
# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Join all three datasets with coordinates
IUCN_with_coords <- CAPTAIN2_IUCN_data %>%
  dplyr::select(PUID, IUCN = Priority) %>%
  left_join(coords, by = "PUID")

EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>%
  dplyr::select(PUID, EDGE2 = Priority) %>%
  left_join(coords, by = "PUID") 

FUSE_with_coords <- CAPTAIN2_FUSE_data %>%
  dplyr::select(PUID, FUSE = Priority) %>%
  left_join(coords, by = "PUID")

# Combine all datasets
all_indices <- coords %>%
  left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by = "PUID")

# Calculate differences
all_indices <- all_indices %>%
  mutate(
    # Replace NA with 0 for calculation purposes
    IUCN = ifelse(is.na(IUCN), 0, IUCN),
    EDGE2 = ifelse(is.na(EDGE2), 0, EDGE2),
    FUSE = ifelse(is.na(FUSE), 0, FUSE),
    
    # Calculate differences
    IUCN_minus_FUSE = IUCN - FUSE,
    IUCN_minus_EDGE2 = IUCN - EDGE2,
    EDGE2_minus_FUSE = EDGE2 - FUSE
  )

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Filter to non-zero differences for each comparison to reduce plot size
# IUCN - FUSE
IUCN_FUSE_diff <- all_indices %>%
  filter(IUCN_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# IUCN - EDGE2
IUCN_EDGE2_diff <- all_indices %>%
  filter(IUCN_minus_EDGE2 != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# EDGE2 - FUSE
EDGE2_FUSE_diff <- all_indices %>%
  filter(EDGE2_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create a diverging color palette for difference maps
# Blue for negative (first index lower), white for zero, red for positive (first index higher)
diff_colors <- c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", 
                "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Display all plots
print(IUCN_FUSE_plot)

Code
print(IUCN_EDGE2_plot)

Code
print(EDGE2_FUSE_plot)

Code
# Save all plots
ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference_2.png"),
  plot = IUCN_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference_2.png"),
  plot = IUCN_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference_2.png"),
  plot = EDGE2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Species level priorities

Code
# Rename column in shark_metrics to match better
shark_metrics <- shark_metrics %>%
  rename(Species = `Species name`)

# Order shark_metrics alphabetically by Species name
shark_metrics <- shark_metrics %>%
  arrange(Species)

# Add an original order ID to protected_fractions to maintain its original order
protected_fractions$original_order <- 1:nrow(protected_fractions)

# Add row number as IDs to both datasets
protected_fractions$protected_ID <- 1:nrow(protected_fractions)
shark_metrics$species_ID <- 1:nrow(shark_metrics)

# Check if the datasets have the same number of rows
if(nrow(protected_fractions) == nrow(shark_metrics)) {
  # Create index mapping
  indices <- data.frame(
    protected_ID = 1:nrow(protected_fractions),
    species_ID = 1:nrow(shark_metrics)
  )
  
  # Join protected_fractions with indices
  protected_with_indices <- protected_fractions %>%
    left_join(indices, by = "protected_ID")
  
  # Join shark_metrics with indices
  shark_with_indices <- shark_metrics %>%
    left_join(indices, by = "species_ID")
  
  # Now join the datasets
  combined_data <- protected_with_indices %>%
    inner_join(
      shark_with_indices,
      by = c("species_ID", "protected_ID"),
      suffix = c("_captain", "_original")
    ) %>%
    arrange(original_order)
  
  cat("Successfully joined datasets with", nrow(combined_data), "species\n")
  
  # Define IUCN categories and order
  iucn_order <- c("LC", "NT", "VU", "EN", "CR")
  
  iucn_colors <- c(
    "LC" = "#50C878",
    "NT" = "#FFFF00",
    "VU" = "#FFA500",
    "EN" = "#FF8C00",
    "CR" = "#FF0000"
  )
  
  # Calculate sample sizes for each group
  iucn_n <- combined_data %>%
    group_by(IUCN_original) %>%
    summarise(n = n()) %>%
    arrange(IUCN_original)
  
  fuse_n <- combined_data %>%
    group_by(FUSE_original) %>%
    summarise(n = n()) %>%
    arrange(FUSE_original)
  
  edge2_n <- combined_data %>%
    group_by(EDGE2_original) %>%
    summarise(n = n()) %>%
    arrange(EDGE2_original)
  
  # Create x-axis labels with sample sizes
  iucn_labels <- paste0(iucn_order, "\n(n=", iucn_n$n, ")")
  fuse_labels <- paste0(1:5, "\n(n=", fuse_n$n, ")")
  edge2_labels <- paste0(1:5, "\n(n=", edge2_n$n, ")")
  
  # 1. IUCN plot
  iucn_plot <- combined_data %>%
    mutate(
      IUCN_status = factor(IUCN_original, levels = 1:5, labels = iucn_order),
      protection_percentage = IUCN_captain * 100
    ) %>%
    ggplot(aes(x = IUCN_status, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "IUCN Red List threat status",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(limits = iucn_order, labels = iucn_labels) +
    scale_fill_manual(values = iucn_colors) +
    scale_color_manual(values = iucn_colors) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 2. FUSE plot
  fuse_plot <- combined_data %>%
    mutate(
      FUSE_category = factor(FUSE_original),
      protection_percentage = FUSE_captain * 100
    ) %>%
    ggplot(aes(x = FUSE_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "FUSE Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = fuse_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 3. EDGE2 plot
  edge2_plot <- combined_data %>%
    mutate(
      EDGE2_category = factor(EDGE2_original),
      protection_percentage = EDGE2_captain * 100
    ) %>%
    ggplot(aes(x = EDGE2_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "EDGE2 Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = edge2_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # Combine the three plots into a grid
  combined_plots <- ggpubr::ggarrange(
    iucn_plot,
    fuse_plot,
    edge2_plot,
    ncol = 3,
    nrow = 1,
    labels = c("(A)", "(B)", "(C)"),
    font.label = list(size = 14,
                      face = "bold",
                      hjust = -2,
                      vjust = 1.5)
  )
  
  print(combined_plots)
  
  # Save combined plot
  ggsave(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"),
         combined_plots, width = 20, height = 6, dpi = 150, bg = "white")
  
  #======================
  # Beta Regression Analysis
  #======================
  
  # Prepare data for beta regression
  model_data <- combined_data %>%
    mutate(
      IUCN_status = IUCN_original,
      FUSE_score = FUSE_original,
      EDGE2_score = EDGE2_original,
      IUCN_protection = IUCN_captain * 100,
      FUSE_protection = FUSE_captain * 100,
      EDGE2_protection = EDGE2_captain * 100,
      # Convert to (0,1) scale and handle boundaries
      n_obs = n(),
      IUCN_protection_beta = (IUCN_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      FUSE_protection_beta = (FUSE_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      EDGE2_protection_beta = (EDGE2_protection/100 * (n_obs - 1) + 0.5) / n_obs
    )
  
  # Fit beta regression models
  iucn_beta <- betareg(IUCN_protection_beta ~ IUCN_status, data = model_data)
  fuse_beta <- betareg(FUSE_protection_beta ~ FUSE_score, data = model_data)
  edge2_beta <- betareg(EDGE2_protection_beta ~ EDGE2_score, data = model_data)
  
  # View summaries
  cat("\n=== IUCN Beta Regression ===\n")
  print(summary(iucn_beta))
  
  cat("\n=== FUSE Beta Regression ===\n")
  print(summary(fuse_beta))
  
  cat("\n=== EDGE2 Beta Regression ===\n")
  print(summary(edge2_beta))
  
  # Get tidy summaries
  iucn_tidy <- tidy(iucn_beta)
  fuse_tidy <- tidy(fuse_beta)
  edge2_tidy <- tidy(edge2_beta)
  
  # Get pseudo R²
  iucn_r2 <- cor(predict(iucn_beta), model_data$IUCN_protection_beta)^2
  fuse_r2 <- cor(predict(fuse_beta), model_data$FUSE_protection_beta)^2
  edge2_r2 <- cor(predict(edge2_beta), model_data$EDGE2_protection_beta)^2
  
  # Calculate marginal effects (average change in percentage per unit increase)
  iucn_margins <- margins(iucn_beta)
  fuse_margins <- margins(fuse_beta)
  edge2_margins <- margins(edge2_beta)
  
  # Extract AME and SE
  iucn_ame_summary <- summary(iucn_margins)
  fuse_ame_summary <- summary(fuse_margins)
  edge2_ame_summary <- summary(edge2_margins)
  
  # Convert to percentage (margins gives proportions)
  iucn_ame <- iucn_ame_summary$AME * 100
  iucn_ame_se <- iucn_ame_summary$SE * 100
  
  fuse_ame <- fuse_ame_summary$AME * 100
  fuse_ame_se <- fuse_ame_summary$SE * 100
  
  edge2_ame <- edge2_ame_summary$AME * 100
  edge2_ame_se <- edge2_ame_summary$SE * 100
  
  cat("\n=== Average Marginal Effects (as percentages) ===\n")
  cat(paste("IUCN AME:", round(iucn_ame, 2), "%, SE:", round(iucn_ame_se, 2), "%\n"))
  cat(paste("FUSE AME:", round(fuse_ame, 2), "%, SE:", round(fuse_ame_se, 2), "%\n"))
  cat(paste("EDGE2 AME:", round(edge2_ame, 2), "%, SE:", round(edge2_ame_se, 2), "%\n"))
  
  # Calculate predicted values and differences
  # IUCN predictions
  iucn_pred <- predict(iucn_beta,
                       newdata = data.frame(IUCN_status = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by IUCN category:\n")
  print(data.frame(Category = 1:5,
                   IUCN_Category = c("LC", "NT", "VU", "EN", "CR"),
                   Predicted_Protection = round(iucn_pred, 2)))
  
  iucn_avg_increase <- mean(diff(iucn_pred))
  cat(paste("Average increase per IUCN category:", round(iucn_avg_increase, 2), "%\n"))
  
  # FUSE predictions
  fuse_pred <- predict(fuse_beta,
                       newdata = data.frame(FUSE_score = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by FUSE score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(fuse_pred, 2)))
  
  fuse_avg_increase <- mean(diff(fuse_pred))
  cat(paste("Average increase per FUSE unit:", round(fuse_avg_increase, 2), "%\n"))
  
  # EDGE2 predictions
  edge2_pred <- predict(edge2_beta,
                        newdata = data.frame(EDGE2_score = 1:5),
                        type = "response") * 100
  
  cat("\nPredicted protection by EDGE2 score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(edge2_pred, 2)))
  
  edge2_avg_increase <- mean(diff(edge2_pred))
  cat(paste("Average increase per EDGE2 unit:", round(edge2_avg_increase, 2), "%\n"))
  
  # Create comprehensive summary table with AME
  comprehensive_summary <- data.frame(
    Approach = c("IUCN", "FUSE", "EDGE2"),
    Coefficient_logit = c(
      iucn_tidy$estimate[2],
      fuse_tidy$estimate[2],
      edge2_tidy$estimate[2]
    ),
    SE = c(
      iucn_tidy$std.error[2],
      fuse_tidy$std.error[2],
      edge2_tidy$std.error[2]
    ),
    z_value = c(
      iucn_tidy$statistic[2],
      fuse_tidy$statistic[2],
      edge2_tidy$statistic[2]
    ),
    P_value = c(
      iucn_tidy$p.value[2],
      fuse_tidy$p.value[2],
      edge2_tidy$p.value[2]
    ),
    AME_percent = c(
      iucn_ame,
      fuse_ame,
      edge2_ame
    ),
    AME_SE = c(
      iucn_ame_se,
      fuse_ame_se,
      edge2_ame_se
    ),
    Pseudo_R2 = c(iucn_r2, fuse_r2, edge2_r2)
  )
  
  cat("\n=== Comprehensive Beta Regression Summary ===\n")
  print(comprehensive_summary)
  
  # Save results
  write.csv(comprehensive_summary,
            here::here("Outputs", "beta_regression_summary.csv"),
            row.names = FALSE)
  
} else {
  cat("ERROR: Datasets have different number of rows.\n")
  cat("Protected fractions:", nrow(protected_fractions), "rows\n")
  cat("Shark metrics:", nrow(shark_metrics), "rows\n")
}
Successfully joined datasets with 1000 species


=== IUCN Beta Regression ===

Call:
betareg(formula = IUCN_protection_beta ~ IUCN_status, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5313 -0.5004 -0.0020  0.5608  2.0202 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.05978    0.07511  -0.796    0.426    
IUCN_status  0.11768    0.02988   3.938 8.21e-05 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.91473    0.03196   28.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 379.5 on 3 Df
Pseudo R-squared: 0.01504
Number of iterations: 7 (BFGS) + 1 (Fisher scoring) 

=== FUSE Beta Regression ===

Call:
betareg(formula = FUSE_protection_beta ~ FUSE_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.3375 -0.4546  0.0051  0.5057  2.3049 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.19664    0.08966  -2.193   0.0283 *
FUSE_score   0.14261    0.06830   2.088   0.0368 *

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  1.04726    0.03673   28.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 211.9 on 3 Df
Pseudo R-squared: 0.004321
Number of iterations: 7 (BFGS) + 2 (Fisher scoring) 

=== EDGE2 Beta Regression ===

Call:
betareg(formula = EDGE2_protection_beta ~ EDGE2_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5490 -0.4954 -0.0052  0.4548  2.1662 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.68964    0.09615  -7.173 7.34e-13 ***
EDGE2_score  0.46406    0.07688   6.036 1.58e-09 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.84174    0.02926   28.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 506.7 on 3 Df
Pseudo R-squared: 0.04761
Number of iterations: 13 (BFGS) + 1 (Fisher scoring) 

=== Average Marginal Effects (as percentages) ===
IUCN AME: 2.9 %, SE: 0.73 %
FUSE AME: 3.56 %, SE: 1.7 %
EDGE2 AME: 11.39 %, SE: 1.85 %

Predicted protection by IUCN category:
  Category IUCN_Category Predicted_Protection
1        1            LC                51.45
2        2            NT                54.38
3        3            VU                57.28
4        4            EN                60.13
5        5            CR                62.92
Average increase per IUCN category: 2.87 %

Predicted protection by FUSE score:
  Score Predicted_Protection
1     1                48.65
2     2                52.21
3     3                55.75
4     4                59.24
5     5                62.63
Average increase per FUSE unit: 3.5 %

Predicted protection by EDGE2 score:
  Score Predicted_Protection
1     1                44.38
2     2                55.93
3     3                66.88
4     4                76.25
5     5                83.63
Average increase per EDGE2 unit: 9.81 %

=== Comprehensive Beta Regression Summary ===
            Approach Coefficient_logit         SE  z_value      P_value
IUCN_status     IUCN         0.1176756 0.02987957 3.938332 8.205012e-05
FUSE_score      FUSE         0.1426108 0.06830262 2.087925 3.680458e-02
EDGE2_score    EDGE2         0.4640609 0.07687772 6.036350 1.576390e-09
            AME_percent    AME_SE  Pseudo_R2
IUCN_status    2.898230 0.7271499 0.01750165
FUSE_score     3.558461 1.6986270 0.00239858
EDGE2_score   11.386873 1.8514148 0.02982912

Manuscript maps

Individual maps

Code
CAPTAIN2_EDGE2_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority EDGE2",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2
CAPTAIN2_EDGE2_msmap

Code
CAPTAIN2_FUSE_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority FUSE",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE
CAPTAIN2_FUSE_msmap

Code
CAPTAIN2_IUCN_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority IUCN",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN
CAPTAIN2_IUCN_msmap

Code
# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid1_maps_continental_2ndrun_new.png"), grid1, width = 8, height = 15, dpi = 300)

Difference maps

Code
# 1. IUCN - FUSE Difference Map
IUCN_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "A")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "B")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid3_maps_continental_2ndrun.png"), grid2, width = 16, height = 15, dpi = 300)

Combined individual and difference maps

Code
library(patchwork)

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Second grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "D")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "E")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "F")

# Create each grid (3 rows, 1 column)
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Combine the two grids side by side (2 columns)
combined_grid <- grid1 | grid2

# Display the combined grid
combined_grid

Code
# If you want to save it
ggsave(
  filename = here::here("Outputs", "Fig.3_combined_priority_maps.png"),
  plot = combined_grid,
  width = 10,    # Adjust width as needed for two columns
  height = 12,   # Adjust height as needed for three rows
  dpi = 300,
  bg = "white"
)

Combined individual and difference maps

Code
library(patchwork)

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Second grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "D")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "E")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "F")

# Create each grid (3 rows, 1 column)
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Combine the two grids side by side (2 columns)
combined_grid <- grid1 | grid2

# Display the combined grid
combined_grid

Code
# If you want to save it
ggsave(
  filename = here::here("outputs", "combined_priority_maps.png"),
  plot = combined_grid,
  width = 10,    # Adjust width as needed for two columns
  height = 12,   # Adjust height as needed for three rows
  dpi = 300,
  bg = "white"
)

Spatial distribution analysis

Code
# Filter for high priority cells (>0.9) and calculate percentages
# Calculate total coastal cells and high-priority cells for each index
# IUCN
CAPTAIN2_IUCN_high <- CAPTAIN2_IUCN_sf %>% 
  filter(Priority > 0.9)
iucn_total_cells <- nrow(CAPTAIN2_IUCN_sf)
iucn_high_cells <- nrow(CAPTAIN2_IUCN_high)
iucn_percent <- round((iucn_high_cells / iucn_total_cells) * 100, 1)

# FUSE
CAPTAIN2_FUSE_high <- CAPTAIN2_FUSE_sf %>% 
  filter(Priority > 0.9)
fuse_total_cells <- nrow(CAPTAIN2_FUSE_sf)
fuse_high_cells <- nrow(CAPTAIN2_FUSE_high)
fuse_percent <- round((fuse_high_cells / fuse_total_cells) * 100, 1)

# EDGE2
CAPTAIN2_EDGE2_high <- CAPTAIN2_EDGE2_sf %>% 
  filter(Priority > 0.9)
edge2_total_cells <- nrow(CAPTAIN2_EDGE2_sf)
edge2_high_cells <- nrow(CAPTAIN2_EDGE2_high)
edge2_percent <- round((edge2_high_cells / edge2_total_cells) * 100, 1)

# Print results
cat("IUCN: ", iucn_high_cells, " cells (", iucn_percent, "%) with priority >0.9\n", sep="")
IUCN: 4734 cells (78.5%) with priority >0.9
Code
cat("FUSE: ", fuse_high_cells, " cells (", fuse_percent, "%) with priority >0.9\n", sep="")
FUSE: 4396 cells (60.5%) with priority >0.9
Code
cat("EDGE2: ", edge2_high_cells, " cells (", edge2_percent, "%) with priority >0.9\n", sep="")
EDGE2: 5328 cells (94.8%) with priority >0.9
Code
#Mean nearest neighbor for high priority grid cells:

# Function to calculate Mean Nearest Neighbor Distance
calculate_mnn <- function(sf_high_priority) {
  # Get centroids and convert to coordinates
  coords <- st_coordinates(st_centroid(sf_high_priority))
  
  # Create a bounding box for the window
  bbox <- st_bbox(sf_high_priority)
  
  # Create observation window
  window <- owin(xrange = c(bbox["xmin"], bbox["xmax"]),
                 yrange = c(bbox["ymin"], bbox["ymax"]))
  
  # Create point pattern object
  ppp_obj <- ppp(coords[,1], coords[,2], window = window)
  
  # Calculate mean nearest neighbor distance
  mnn <- mean(nndist(ppp_obj))
  
  # Calculate standard deviation for context
  mnn_sd <- sd(nndist(ppp_obj))
  
  return(list(mean = mnn, sd = mnn_sd))
}

# Calculate for each index
iucn_mnn <- calculate_mnn(CAPTAIN2_IUCN_high)
fuse_mnn <- calculate_mnn(CAPTAIN2_FUSE_high)
edge2_mnn <- calculate_mnn(CAPTAIN2_EDGE2_high)

# Print results
cat("\nMean Nearest Neighbor Distance (units depend on your projection):\n")

Mean Nearest Neighbor Distance (units depend on your projection):
Code
cat("IUCN:  Mean = ", round(iucn_mnn$mean, 2), ", SD = ", round(iucn_mnn$sd, 2), "\n", sep="")
IUCN:  Mean = 51143.37, SD = 17182.87
Code
cat("FUSE:  Mean = ", round(fuse_mnn$mean, 2), ", SD = ", round(fuse_mnn$sd, 2), "\n", sep="")
FUSE:  Mean = 49594.94, SD = 26765.94
Code
cat("EDGE2: Mean = ", round(edge2_mnn$mean, 2), ", SD = ", round(edge2_mnn$sd, 2), "\n", sep="")
EDGE2: Mean = 48211.63, SD = 8975.76
Code
# Create summary data frame
mnn_summary <- data.frame(
  Index = c("IUCN", "FUSE", "EDGE2"),
  N_cells = c(iucn_high_cells, fuse_high_cells, edge2_high_cells),
  Percent = c(iucn_percent, fuse_percent, edge2_percent),
  MNN_mean = c(iucn_mnn$mean, fuse_mnn$mean, edge2_mnn$mean),
  MNN_sd = c(iucn_mnn$sd, fuse_mnn$sd, edge2_mnn$sd)
)

print(mnn_summary)
  Index N_cells Percent MNN_mean    MNN_sd
1  IUCN    4734    78.5 51143.37 17182.874
2  FUSE    4396    60.5 49594.94 26765.936
3 EDGE2    5328    94.8 48211.63  8975.762

Congruence maps

Code
# Method 1: If you have separate dataframes for each index
# Assuming your data has columns: PUID, Priority, and geometry

# Fxirst, identify high priority areas (>0.9) for each index
high_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority > 0.9, ]
high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority > 0.9, ]
high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority > 0.9, ]

# Find congruent areas (present in all three)
# Method using PUID (assuming you have PUID column)
congruent_PUIDs <- intersect(intersect(high_priority_EDGE2$PUID, 
                                      high_priority_FUSE$PUID), 
                            high_priority_IUCN$PUID)

# Extract congruent areas
congruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]

# Print summary
cat("High priority areas (>0.9):\n")
High priority areas (>0.9):
Code
cat("EDGE2:", nrow(high_priority_EDGE2), "areas\n")
EDGE2: 5328 areas
Code
cat("FUSE:", nrow(high_priority_FUSE), "areas\n") 
FUSE: 4396 areas
Code
cat("IUCN:", nrow(high_priority_IUCN), "areas\n")
IUCN: 4734 areas
Code
cat("Congruent areas:", nrow(congruent_areas), "areas\n")
Congruent areas: 1459 areas
Code
cat("Percentage of overlap:", round(nrow(congruent_areas)/min(nrow(high_priority_EDGE2), 
                                                            nrow(high_priority_FUSE), 
                                                            nrow(high_priority_IUCN))*100, 2), "%\n")
Percentage of overlap: 33.19 %
Code
# Method 2: If you need to merge dataframes first
# Create a combined dataframe with all three priority scores
# First, extract just the data (without geometry) from the other sf objects
FUSE_data <- st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])
IUCN_data <- st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])

# Merge with the EDGE2 sf object (keeping geometry)
combined_priorities <- merge(CAPTAIN2_EDGE2_sf, FUSE_data, 
                            by = "PUID", suffixes = c("_EDGE2", "_FUSE"))
combined_priorities <- merge(combined_priorities, IUCN_data, 
                            by = "PUID")
names(combined_priorities)[names(combined_priorities) == "Priority"] <- "Priority_IUCN"

# Identify congruent high priority areas
congruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 > 0.9 & 
                                         combined_priorities$Priority_FUSE > 0.9 & 
                                         combined_priorities$Priority_IUCN > 0.9, ]

# Create a map showing only congruent areas
congruent_map <- ggplot() +
  geom_sf(data = congruent_areas, aes(fill = "Congruent High Priority"), 
          color = "red", size = 0.8, alpha = 0.8) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_fill_manual(values = c("Congruent high priority" = "red"),
                    name = "Priority areas") +
  labs(title = "Congruent high priority areas",
       subtitle = "Areas with priority > 0.9 in IUCN, FUSE and EDGE2 dimensions",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the map
print(congruent_map)

Code
ggsave(here::here("Outputs", "congruent_high_prioritiy_areas_09.png"), 
       congruent_map, 
       width = 10, height = 8, dpi = 300, bg = "white")

# Create detailed congruence categories
# Define logical conditions for each index
EDGE2_high <- combined_priorities$Priority_EDGE2 > 0.9
FUSE_high <- combined_priorities$Priority_FUSE > 0.9
IUCN_high <- combined_priorities$Priority_IUCN > 0.9

# Create detailed congruence categories
combined_priorities$congruence_category <- case_when(
  EDGE2_high & FUSE_high & IUCN_high ~ "All three indices",
  !EDGE2_high & FUSE_high & IUCN_high ~ "IUCN + FUSE",      # Changed from "FUSE + IUCN"
  EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN + EDGE2",     # Changed from "EDGE2 + IUCN"
  EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE + EDGE2",     # Changed from "EDGE2 + FUSE"
  !EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN only",
  !EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE only",
  EDGE2_high & !FUSE_high & !IUCN_high ~ "EDGE2 only",
  TRUE ~ "No high priority"
)

# Convert to factor with desired order
combined_priorities$congruence_category <- factor(
  combined_priorities$congruence_category,
  levels = c("All three indices",
             "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed order
             "IUCN only", "FUSE only", "EDGE2 only")
)

# Filter data to only include high priority areas
high_priority_data <- combined_priorities[combined_priorities$congruence_category %in%
                                          c("All three indices", 
                                            "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed names
                                            "IUCN only", "FUSE only", "EDGE2 only"), ]

# Debug: Check the data
cat("Data check:\n")
Data check:
Code
cat("Total high priority areas:", nrow(high_priority_data), "\n")
Total high priority areas: 3139 
Code
print(table(high_priority_data$congruence_category))

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
             1459                75              1130               275 
        IUCN only         FUSE only        EDGE2 only 
               29                17               154 
Code
# Check if the factor levels are properly set
cat("\nFactor levels:\n")

Factor levels:
Code
print(levels(high_priority_data$congruence_category))
[1] "All three indices" "IUCN + FUSE"       "IUCN + EDGE2"     
[4] "FUSE + EDGE2"      "IUCN only"         "FUSE only"        
[7] "EDGE2 only"       
Code
# Check for any issues with the data
cat("\nData structure check:\n")

Data structure check:
Code
print(str(high_priority_data$congruence_category))
 Factor w/ 7 levels "All three indices",..: 6 6 6 6 6 6 6 6 4 4 ...
NULL
Code
# Test a simple plot first
cat("\nTesting simple plot...\n")

Testing simple plot...
Code
test_plot <- ggplot() +
  geom_sf(data = high_priority_data, aes(fill = congruence_category)) +
  scale_fill_manual(values = c("All three indices" = "#8B0000",
                              "EDGE2 + FUSE" = "#FF4500",
                              "EDGE2 + IUCN" = "#FF6347",
                              "FUSE + IUCN" = "#FFA500",
                              "EDGE2 only" = "#4169E1",
                              "FUSE only" = "#32CD32",
                              "IUCN only" = "#9370DB")) +
  theme_void()

# Try to print the simple test plot
tryCatch({
  print(test_plot)
  cat("Simple plot works!\n")
}, error = function(e) {
  cat("Simple plot failed with error:", e$message, "\n")
})

Simple plot works!
Code
# If simple plot fails, let's check the data more thoroughly
if (nrow(high_priority_data) == 0) {
  cat("No high priority data found! Checking original data...\n")
  cat("EDGE2 > 0.9:", sum(combined_priorities$Priority_EDGE2 > 0.9, na.rm = TRUE), "\n")
  cat("FUSE > 0.9:", sum(combined_priorities$Priority_FUSE > 0.9, na.rm = TRUE), "\n") 
  cat("IUCN > 0.9:", sum(combined_priorities$Priority_IUCN > 0.9, na.rm = TRUE), "\n")
}

# Extract coordinates for the congruence plot
if (nrow(high_priority_data) > 0) {
  coords <- st_coordinates(st_centroid(high_priority_data))
  plot_data <- data.frame(
    x = coords[,1],
    y = coords[,2], 
    category = high_priority_data$congruence_category
  )
  
  # Remove any rows with NA category
  plot_data <- plot_data[!is.na(plot_data$category), ]
  
  cat("Creating enhanced congruence plot...\n")
  cat("Plot data dimensions:", nrow(plot_data), "points\n")
  
  # Create the enhanced congruence plot
  congruence_map <- ggplot(plot_data, aes(x = x, y = y, color = category)) +
    geom_point(size = 1.2, alpha = 0.85, stroke = 0) +  # Larger points, no stroke to reduce overlap
    scale_color_manual(
      values = c("All three indices" = "#8B0000",      # Dark red
                "EDGE2 + FUSE" = "#FF8C00",            # Dark orange  
                "EDGE2 + IUCN" = "#DC143C",            # Crimson
                "FUSE + IUCN" = "#FFD700",             # Gold
                "EDGE2 only" = "#4169E1",              # Royal blue
                "FUSE only" = "#32CD32",               # Lime green
                "IUCN only" = "#9370DB"),              # Medium purple
      name = "High Priority\nCongruence",
      guide = guide_legend(
        override.aes = list(size = 4, alpha = 1),  # Larger, more opaque legend points
        title.position = "top",
        title.hjust = 0.5,
        ncol = 4  # Two columns for legend to save space
      )) +
    labs(
      title = "Global Conservation Priority Congruence Analysis",
      subtitle = "High priority areas (>0.9) showing agreement patterns across EDGE2, FUSE, and IUCN indices",
      x = "Longitude", 
      y = "Latitude",
      caption = paste("Total high priority areas:", nrow(plot_data))
    ) +
    theme_minimal() +
    theme(
      # Plot aesthetics
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "#f8f9fa", color = NA),
      panel.grid.major = element_line(color = "white", size = 0.5, linetype = "solid"),
      panel.grid.minor = element_line(color = "white", size = 0.25, linetype = "solid"),
      
      # Text styling
      plot.title = element_text(size = 16, hjust = 0.5, face = "bold", 
                               margin = margin(b = 10)),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray30",
                                  margin = margin(b = 20)),
      plot.caption = element_text(size = 10, color = "gray50", hjust = 1),
      
      # Axis styling
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(size = 10, color = "gray30"),
      axis.ticks = element_line(color = "gray50", size = 0.5),
      
      # Legend styling
      legend.position = "bottom",
      legend.background = element_rect(fill = "white", color = "gray80", size = 0.5),
      legend.margin = margin(t = 15, r = 10, b = 10, l = 10),
      legend.title = element_text(size = 12, face = "bold"),
      legend.text = element_text(size = 10),
      legend.key.size = unit(0.8, "cm"),
      
      # Panel border
      panel.border = element_rect(color = "gray70", fill = NA, size = 0.5)
    ) +
    # Add coordinate system and aspect ratio
    coord_fixed(ratio = 1.3) +  # Adjust ratio for better map appearance
    # Add subtle borders around plot area
    scale_x_continuous(expand = c(0.02, 0.02)) +
    scale_y_continuous(expand = c(0.02, 0.02))
  
  # Print the enhanced plot
  print(congruence_map)
  cat("Enhanced congruence plot created successfully!\n")
  
  # Print summary by category
  cat("\nBreakdown by congruence category:\n")
  category_counts <- table(plot_data$category)
  category_percentages <- round(prop.table(category_counts) * 100, 1)
  
  for(i in 1:length(category_counts)) {
    cat(sprintf("%-20s: %4d points (%4.2f%%)\n", 
                names(category_counts)[i], 
                category_counts[i], 
                category_percentages[i]))
  }
  
} else {
  cat("No valid high priority data to plot\n")
}
Creating enhanced congruence plot...
Plot data dimensions: 3139 points

Enhanced congruence plot created successfully!

Breakdown by congruence category:
All three indices   : 1459 points (46.50%)
IUCN + FUSE         :   75 points (2.40%)
IUCN + EDGE2        : 1130 points (36.00%)
FUSE + EDGE2        :  275 points (8.80%)
IUCN only           :   29 points (0.90%)
FUSE only           :   17 points (0.50%)
EDGE2 only          :  154 points (4.90%)
Code
# Print summary of congruence patterns
cat("\nCongruence Pattern Summary:\n")

Congruence Pattern Summary:
Code
congruence_summary <- table(combined_priorities$congruence_category, useNA = "ifany")
print(congruence_summary)

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
             1459                75              1130               275 
        IUCN only         FUSE only        EDGE2 only              <NA> 
               29                17               154                 6 
Code
# Calculate percentages
total_high_priority <- sum(congruence_summary[names(congruence_summary) != "No high priority"], na.rm = TRUE)
congruence_percentages <- round(congruence_summary / total_high_priority * 100, 2)
cat("\nPercentages of high priority areas:\n")

Percentages of high priority areas:
Code
print(congruence_percentages[names(congruence_percentages) != "No high priority"])

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
            46.48              2.39             36.00              8.76 
        IUCN only         FUSE only        EDGE2 only              <NA> 
             0.92              0.54              4.91                   
Code
# ms style map -----
# Step 1: Check original data
print("=== STEP 1: CHECKING ORIGINAL DATA ===")
[1] "=== STEP 1: CHECKING ORIGINAL DATA ==="
Code
print(paste("Number of rows in plot_data:", nrow(plot_data)))
[1] "Number of rows in plot_data: 3139"
Code
print("First few rows of plot_data:")
[1] "First few rows of plot_data:"
Code
print(head(plot_data))
          x       y  category
1 -767764.9 5762063 FUSE only
2 -727356.2 5762063 FUSE only
3 -772256.2 5713698 FUSE only
4 -731611.1 5713698 FUSE only
5 -776677.2 5665099 FUSE only
6 -735799.5 5665099 FUSE only
Code
print("Summary of coordinates:")
[1] "Summary of coordinates:"
Code
print(summary(plot_data[c("x", "y")]))
       x                   y           
 Min.   :-10914637   Min.   :-5217532  
 1st Qu.: -3784462   1st Qu.:-1511094  
 Median :  8930102   Median :  848691  
 Mean   :  4702046   Mean   :  599936  
 3rd Qu.: 12115021   3rd Qu.: 2810253  
 Max.   : 16383151   Max.   : 5762063  
Code
print("Unique categories:")
[1] "Unique categories:"
Code
print(unique(plot_data$category))
[1] FUSE only         FUSE + EDGE2      All three indices EDGE2 only       
[5] IUCN + EDGE2      IUCN + FUSE       IUCN only        
7 Levels: All three indices IUCN + FUSE IUCN + EDGE2 ... EDGE2 only
Code
# Step 2: Transform to sf object
print("=== STEP 2: CREATING SF OBJECT ===")
[1] "=== STEP 2: CREATING SF OBJECT ==="
Code
congruence_sf <- st_as_sf(
  plot_data, 
  coords = c("x", "y"), 
  crs = mcbryde_thomas_2  # Assuming WGS84
)
print(paste("Number of sf features created:", nrow(congruence_sf)))
[1] "Number of sf features created: 3139"
Code
print("Original CRS:")
[1] "Original CRS:"
Code
print(st_crs(congruence_sf))
Coordinate Reference System:
  User input: +proj=mbt_s 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ mbt_s"]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Code
print("Bounding box before projection:")
[1] "Bounding box before projection:"
Code
print(st_bbox(congruence_sf))
     xmin      ymin      xmax      ymax 
-10914637  -5217532  16383151   5762063 
Code
# Step 3: Project to McBryde-Thomas 2
print("=== STEP 3: PROJECTING TO MCBRYDE-THOMAS 2 ===")
[1] "=== STEP 3: PROJECTING TO MCBRYDE-THOMAS 2 ==="
Code
congruence_sf <- st_transform(congruence_sf, crs = mcbryde_thomas_2)
print("CRS after projection:")
[1] "CRS after projection:"
Code
print(st_crs(congruence_sf))
Coordinate Reference System:
  User input: +proj=mbt_s 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ mbt_s"]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Code
print("Bounding box after projection:")
[1] "Bounding box after projection:"
Code
print(st_bbox(congruence_sf))
     xmin      ymin      xmax      ymax 
-10914637  -5217532  16383151   5762063 
Code
# Step 4: Check world projection
print("=== STEP 4: CHECKING WORLD PROJECTION ===")
[1] "=== STEP 4: CHECKING WORLD PROJECTION ==="
Code
world_projected_congruence <- st_transform(world, crs = mcbryde_thomas_2)
print("World bounding box:")
[1] "World bounding box:"
Code
print(st_bbox(world_projected_congruence))
     xmin      ymin      xmax      ymax 
-18030969  -8669729  18199360   8324207 
Code
# Step 5: Create globe border
print("=== STEP 5: CREATING GLOBE BORDER ===")
[1] "=== STEP 5: CREATING GLOBE BORDER ==="
Code
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))
globe_border_congruence <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)
print("Globe border bounding box:")
[1] "Globe border bounding box:"
Code
print(st_bbox(globe_border_congruence))
     xmin      ymin      xmax      ymax 
-18373133  -8669783  18373133   8669783 
Code
# Step 6: Check intersection
print("=== STEP 6: CHECKING INTERSECTION ===")
[1] "=== STEP 6: CHECKING INTERSECTION ==="
Code
points_in_globe <- st_intersects(congruence_sf, globe_border_congruence, sparse = FALSE)
print(paste("Points within globe boundary:", sum(points_in_globe)))
[1] "Points within globe boundary: 3139"
Code
print(paste("Percentage of points in globe:", round(sum(points_in_globe)/nrow(congruence_sf)*100, 2), "%"))
[1] "Percentage of points in globe: 100 %"
Code
# Step 7: Extract coordinates for verification
print("=== STEP 7: EXTRACTING COORDINATES ===")
[1] "=== STEP 7: EXTRACTING COORDINATES ==="
Code
coords <- st_coordinates(congruence_sf)
print("Projected coordinate summary:")
[1] "Projected coordinate summary:"
Code
print(summary(coords))
       X                   Y           
 Min.   :-10914637   Min.   :-5217532  
 1st Qu.: -3784462   1st Qu.:-1511094  
 Median :  8930102   Median :  848691  
 Mean   :  4702046   Mean   :  599936  
 3rd Qu.: 12115021   3rd Qu.: 2810253  
 Max.   : 16383151   Max.   : 5762063  
Code
print("First few projected coordinates:")
[1] "First few projected coordinates:"
Code
print(head(coords))
             X       Y
[1,] -767764.9 5762063
[2,] -727356.2 5762063
[3,] -772256.2 5713698
[4,] -731611.1 5713698
[5,] -776677.2 5665099
[6,] -735799.5 5665099
Code
# Step 8: Create theme
my_theme_congruence <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Step 9: Create plot with debugging
print("=== STEP 9: CREATING PLOT ===")
[1] "=== STEP 9: CREATING PLOT ==="
Code
congruence_map <- ggplot() +
  geom_sf(data = congruence_sf, aes(color = category), size = 1.2, alpha = 0.85) +
  geom_sf(data = world_projected_congruence, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_congruence, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_manual(
  values = c("All three indices" = "#8B0000",
            "IUCN + FUSE" = "#FFD700",           # Changed from "FUSE + IUCN"
            "IUCN + EDGE2" = "#DC143C",          # Changed from "EDGE2 + IUCN"
            "FUSE + EDGE2" = "#FF8C00",          # Changed from "EDGE2 + FUSE"
            "IUCN only" = "#9370DB",
            "FUSE only" = "#32CD32",
            "EDGE2 only" = "#4169E1"),
  breaks = c("All three indices",
             "IUCN + FUSE", 
             "IUCN + EDGE2", 
             "FUSE + EDGE2",
             "IUCN only", 
             "FUSE only", 
             "EDGE2 only"),              # Finally: FUSE + EDGE2
    name = "High priority congruence",
    guide = guide_legend(
      override.aes = list(size = 4, alpha = 1),
      title.position = "top", 
      title.hjust = 0.5,
      ncol = 4
    )) +
  labs(
    title = "Global conservation priority congruence analysis",
    subtitle = "High priority areas (>0.9) showing agreement patterns across IUCN, FUSE and EDGE2 prioritisations",
    x = NULL, 
    y = NULL
  ) +
  my_theme_congruence

ggsave(here::here("outputs", "congruent_high_prioritiy_areas_09_all_indices.png"), 
       congruence_map, 
       width = 10, height = 8, dpi = 300, bg = "white")

print("=== STEP 10: DISPLAYING PLOT ===")
[1] "=== STEP 10: DISPLAYING PLOT ==="
Code
print("Plot created successfully")
[1] "Plot created successfully"
Code
# Alternative test plot with just points to see if they're visible
print("=== CREATING TEST PLOT WITH JUST POINTS ===")
[1] "=== CREATING TEST PLOT WITH JUST POINTS ==="
Code
test_plot <- ggplot() +
  geom_sf(data = congruence_sf, aes(color = category), size = 2, alpha = 1) +
  scale_color_manual(
    values = c("All three indices" = "#8B0000",      
              "EDGE2 + FUSE" = "#FF8C00",            
              "EDGE2 + IUCN" = "#DC143C",            
              "FUSE + IUCN" = "#FFD700",             
              "EDGE2 only" = "#4169E1",              
              "FUSE only" = "#32CD32",               
              "IUCN only" = "#9370DB")) +
  labs(title = "Test Plot - Points Only")

print("Test plot created")
[1] "Test plot created"
Code
# Display both plots
print(congruence_map)

Code
print(test_plot)

Congruence tests

Code
# Statistical congruence analysis following the methods section
# First, extract priority values for the comparison
priority_threshold <- 0.9

# Create binary matrices for high-priority areas (>0.9)
fuse_high <- CAPTAIN2_FUSE_data_nonzero$Priority > priority_threshold
edge2_high <- CAPTAIN2_EDGE2_data_nonzero$Priority > priority_threshold
iucn_high <- CAPTAIN2_IUCN_data_nonzero$Priority > priority_threshold

# Calculate observed overlaps
observed_fuse_edge2 <- sum(fuse_high & edge2_high)
observed_fuse_iucn <- sum(fuse_high & iucn_high)
observed_edge2_iucn <- sum(edge2_high & iucn_high)
observed_all_three <- sum(fuse_high & edge2_high & iucn_high)

# Calculate minimum high-priority cells for each comparison
min_fuse_edge2 <- min(sum(fuse_high), sum(edge2_high))
min_fuse_iucn <- min(sum(fuse_high), sum(iucn_high))
min_edge2_iucn <- min(sum(edge2_high), sum(iucn_high))

# Calculate percentage overlaps (conservative estimate)
overlap_fuse_edge2 <- (observed_fuse_edge2 / min_fuse_edge2) * 100
overlap_fuse_iucn <- (observed_fuse_iucn / min_fuse_iucn) * 100
overlap_edge2_iucn <- (observed_edge2_iucn / min_edge2_iucn) * 100

# Randomization test function
randomization_test <- function(index1_high, index2_high, n_permutations = 999) {
  observed_overlap <- sum(index1_high & index2_high)
  min_cells <- min(sum(index1_high), sum(index2_high))
  observed_percentage <- (observed_overlap / min_cells) * 100
  
  # Generate random overlaps
  random_overlaps <- replicate(n_permutations, {
    # Randomly shuffle the spatial distribution of index1
    shuffled_index1 <- sample(index1_high)
    random_overlap <- sum(shuffled_index1 & index2_high)
    (random_overlap / min_cells) * 100
  })
  
  # Calculate p-value
  if (observed_percentage > mean(random_overlaps)) {
    p_value <- sum(random_overlaps >= observed_percentage) / n_permutations
  } else {
    p_value <- sum(random_overlaps <= observed_percentage) / n_permutations
  }
  
  return(list(
    observed_percentage = observed_percentage,
    expected_percentage = mean(random_overlaps),
    p_value = p_value,
    random_overlaps = random_overlaps
  ))
}

# Perform randomization tests
print("=== STATISTICAL CONGRUENCE ANALYSIS ===")
[1] "=== STATISTICAL CONGRUENCE ANALYSIS ==="
Code
cat("Threshold for high-priority areas:", priority_threshold, "\n\n")
Threshold for high-priority areas: 0.9 
Code
# FUSE vs EDGE2
test_fuse_edge2 <- randomization_test(fuse_high, edge2_high)

# FUSE vs IUCN
test_fuse_iucn <- randomization_test(fuse_high, iucn_high)

# EDGE2 vs IUCN
test_edge2_iucn <- randomization_test(edge2_high, iucn_high)

# Create a summary table
congruence_results <- data.frame(
  Comparison = c("FUSE vs EDGE2", "FUSE vs IUCN", "EDGE2 vs IUCN"),
  Observed_Overlap_Percent = c(test_fuse_edge2$observed_percentage,
                               test_fuse_iucn$observed_percentage,
                               test_edge2_iucn$observed_percentage),
  Expected_Overlap_Percent = c(test_fuse_edge2$expected_percentage,
                               test_fuse_iucn$expected_percentage,
                               test_edge2_iucn$expected_percentage),
  P_Value = c(test_fuse_edge2$p_value,
              test_fuse_iucn$p_value,
              test_edge2_iucn$p_value),
  Significance = c(
    ifelse(test_fuse_edge2$p_value < 0.001, "***",
           ifelse(test_fuse_edge2$p_value < 0.01, "**",
                  ifelse(test_fuse_edge2$p_value < 0.05, "*", "ns"))),
    ifelse(test_fuse_iucn$p_value < 0.001, "***",
           ifelse(test_fuse_iucn$p_value < 0.01, "**",
                  ifelse(test_fuse_iucn$p_value < 0.05, "*", "ns"))),
    ifelse(test_edge2_iucn$p_value < 0.001, "***",
           ifelse(test_edge2_iucn$p_value < 0.01, "**",
                  ifelse(test_edge2_iucn$p_value < 0.05, "*", "ns")))
  )
)

print(congruence_results)
     Comparison Observed_Overlap_Percent Expected_Overlap_Percent     P_Value
1 FUSE vs EDGE2                 94.56324                 94.90466 0.062062062
2  FUSE vs IUCN                 74.34031                 75.63198 0.002002002
3 EDGE2 vs IUCN                 95.18378                 94.83053 0.009009009
  Significance
1           ns
2           **
3           **

Relationship with fishing effort : bivariate maps

Code
# Load your CAPTAIN data
IUCN_data <- CAPTAIN2_IUCN_data_nonzero  # Assuming this exists from your previous code
FUSE_data <- CAPTAIN2_FUSE_data_nonzero  # Assuming this exists from your previous code
EDGE2_data <- CAPTAIN2_EDGE2_data_nonzero # Assuming this exists from your previous code

# Process and prepare the data function
process_bivariate_data <- function(priority_data, fishing_data, min_priority = 0) {
  # Filter for cells with priority > min_priority
  priority_data_filtered <- priority_data %>%
    filter(Priority > min_priority)
  
  # Prepare fishing data
  fishing_data <- fishing_data %>%
    rename(Longitude = lon_05deg, Latitude = lat_05deg, FishingHours = predicted_fishing_hours)
  
  # Join datasets
  combined_data <- priority_data_filtered %>%
    left_join(fishing_data, by = c("Longitude", "Latitude"))
  
  # Handle NAs in fishing hours (replace with 0)
  combined_data$FishingHours[is.na(combined_data$FishingHours)] <- 0
  
  # Normalize priorities to 0-1 range and log transform fishing hours
  combined_data <- combined_data %>%
    mutate(
      Priority_Norm = Priority, # Assuming already in 0-1 range
      # Log transform fishing hours to better handle skewed distribution
      FishingHours_Log = log1p(FishingHours),
      # Normalize log fishing hours
      FishingHours_Norm = (FishingHours_Log - min(FishingHours_Log, na.rm = TRUE)) /
                         (max(FishingHours_Log, na.rm = TRUE) - min(FishingHours_Log, na.rm = TRUE))
    )
  
  # Set projection
  mcbryde_thomas_2 <- "+proj=mbt_s"
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process the three datasets with fishing data (only cells with Priority > 0)
iucn_fishing_sf <- process_bivariate_data(IUCN_data, aggregated_data, min_priority = 0)
fuse_fishing_sf <- process_bivariate_data(FUSE_data, aggregated_data, min_priority = 0)
edge2_fishing_sf <- process_bivariate_data(EDGE2_data, aggregated_data, min_priority = 0)

# Create color palette - using the PurpleOr scheme
create_bivariate_palette <- function() {
  map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
  map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
  map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
  map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
  map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
  map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
  map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
  map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
  map_pal_mtx[1, 1] <- '#ffffee'
  map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))
  return(map_pal)
}

map_pal <- create_bivariate_palette()

# Color mapping function
get_color <- function(priority, fishing) {
  priority_class <- cut(priority, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fishing_class <- cut(fishing, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fishing_class)-1)*4 + as.numeric(priority_class)])
}

# Apply colors to the datasets
iucn_fishing_sf$bivariate_color <- mapply(get_color,
                                          iucn_fishing_sf$Priority_Norm,
                                          iucn_fishing_sf$FishingHours_Norm)

fuse_fishing_sf$bivariate_color <- mapply(get_color,
                                          fuse_fishing_sf$Priority_Norm,
                                          fuse_fishing_sf$FishingHours_Norm)

edge2_fishing_sf$bivariate_color <- mapply(get_color,
                                          edge2_fishing_sf$Priority_Norm,
                                          edge2_fishing_sf$FishingHours_Norm)

# Get world map and project
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create globe border
lon_points <- rep(c(-180, 180), each = 100)
lat_points <- c(seq(-90, 90, length.out = 100), seq(90, -90, length.out = 100))
globe_outline <- data.frame(lon = lon_points, lat = lat_points) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_combine() %>%
  st_cast("POLYGON") %>%
  st_transform(crs = mcbryde_thomas_2)

# Function to create individual bivariate maps with labels and titles
create_bivariate_map <- function(data_sf, label, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = bivariate_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_outline, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL, title = title) +  # Added title
    theme(panel.grid = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5),  # Title styling
          plot.margin = margin(5, 5, 5, 5, "mm")) +
    # Add smaller label in top-left corner
    annotate("text", x = -Inf, y = Inf, label = label,
             hjust = -0.2, vjust = 1.2, size = 4, fontface = "bold")  # Reduced from 8 to 4
}

# Create the three maps with labels and titles
iucn_plot <- create_bivariate_map(iucn_fishing_sf, "(A)", "IUCN")
fuse_plot <- create_bivariate_map(fuse_fishing_sf, "(B)", "FUSE")  
edge2_plot <- create_bivariate_map(edge2_fishing_sf, "(C)", "EDGE2")

# Create a larger legend for better visibility
larger_legend <- bi_legend(pal = map_pal, dim = 4,
                          xlab = 'Conservation\npriority',
                          ylab = 'Fishing\neffort',
                          size = 4) +  # Increased from 3 to 4 for bigger legend
  theme(
    axis.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    axis.text = element_blank(),
    legend.text = element_text(size = 8),  # Reduced from 10 to 8
    legend.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    plot.margin = margin(10, 10, 10, 10, "mm")
  )

# Create layout matrix: 3 rows, 2 columns
# Column 1: maps (wider)
# Column 2: empty, legend (middle row), empty
layout_matrix <- rbind(
  c(1, NA),      # Row 1: Map A, empty
  c(2, 4),       # Row 2: Map B, legend
  c(3, NA)       # Row 3: Map C, empty
)

# Arrange all plots together
combined_plot <- grid.arrange(
  iucn_plot,
  fuse_plot,
  edge2_plot,
  larger_legend,
  layout_matrix = layout_matrix,
  heights = c(1, 1, 1),      # Equal height for all three map rows
  widths = c(3, 1)           # Maps take 3/4 width, legend takes 1/4
)

Code
#To display the plot in the quarto file 
grid.arrange(
  iucn_plot,
  fuse_plot,
  edge2_plot,
  larger_legend,
  layout_matrix = layout_matrix,
  heights = c(1, 1, 1),      # Equal height for all three map rows
  widths = c(3, 1)           # Maps take 3/4 width, legend takes 1/4
)

Code
# Save with dimensions better suited for A4 portrait
ggsave(here::here("Outputs", "Fig.4_All_indices_fishing_bivariate_maps_vertical.png"),
       combined_plot, width = 8.3, height = 11, dpi = 300, bg = "white")  # A4 portrait dimensions